# This code replicates the analyses of LAPOP data, producing Figures 7 and A18 and Table A11.

library(dplyr)
library(magrittr)
library(xtable)
library(ggplot2)
library(patchwork)

set.seed(416) ## use this to reproduce bootstrapped estimates

load("data/lapop_extract.Rdata")

# select relevant survey waves in which trust/attendance questions were asked
waves <- dplyr::group_by(lapop_sset, pais_num, wave) %>%
  dplyr::summarize(
    mun_trust_ctr = mean(is.na(mun_trust)),
    comm_trust_ctr = mean(is.na(comm_trust))
  ) %>%
  filter(wave > 2004 & pais_num < 27)

id_list <- filter(lapop_sset, year > 2004 & pais_num < 27) %>% 
  dplyr::group_by(ID) %>%
  dplyr::summarize(
    mun_trust = mean(is.na(mun_trust)),
    mtg = mean(is.na(mtg)),
    comm_trust = mean(is.na(comm_trust)),
    comm_assoc = mean(is.na(comm_assoc))
  ) 

mun_list = id_list$ID[id_list$mtg < 1 & id_list$mun_trust < 1]
comm_list = id_list$ID[id_list$comm_trust < 1 & id_list$comm_assoc < 1]

# splits sample by country-wave (ID)
ind_id <- split(1:nrow(lapop_sset), lapop_sset$ID)

key <- dplyr::group_by(lapop_sset, pais_num, pais, year) %>%
  summarize() %>%
  filter(year > 2004 & pais_num < 27) %>%
  mutate(ID = paste0(pais_num, "_",year)) 


# (block) bootstrap function resamples (by survey wave) and calculates 
bootstrap <- function(){
  # generates a bootstrapped sample
  samp_id <- lapply(ind_id, FUN = function(x){sample(x, size = length(x), replace = T)})
  samp <- lapop_sset[unlist(samp_id),]
  print("sample created")
  # estimates Pr(Attended|Trust = highest) - Pr(Attended|Trust = lowest) for city council mtgs (linearly)
  mun <- sapply(mun_list, FUN = function(x){
    summary(lm(mtg ~ mun_trust, data = samp, subset = ID == x))$coef[2,1]
  })
  # estimates Pr(Attended) for city council meetings
  mun_level <- tapply(samp$mtg, INDEX = samp$ID, FUN = mean, na.rm = T)
  mun_level <- mun_level[names(mun_level) %in% mun_list]
  # estimates E[Trust in local government]
  mun_trust <- tapply(samp$mun_trust, INDEX = samp$ID, FUN = mean, na.rm = T)
  mun_trust <- mun_trust[names(mun_trust) %in% mun_list]
   # estimates Pr(Attended|Trust = highest) - Pr(Attended|Trust = lowest) for community association mtgs (linearly)
  comm <- sapply(comm_list, FUN = function(x){
    summary(lm(comm_assoc ~ comm_trust, data = samp, subset = ID == x))$coef[2,1]
  })
  # estimates Pr(Attended) for community association meetings
  comm_level <- tapply(samp$comm_assoc, INDEX = samp$ID, FUN = mean, na.rm = T)
  comm_level <- comm_level[names(comm_level) %in% comm_list]
  # estimates E[Trust in community]
  comm_trust <- tapply(samp$comm_trust, INDEX = samp$ID, FUN = mean, na.rm = T)
  comm_trust <- comm_trust[names(comm_trust) %in% comm_list]
  
  return(list(mun, mun_level, mun_trust, comm, comm_level, comm_trust))
  
}

# this function executes the bootstrapping and saves estimates
iterate <- function(n = 500){
  reps <- replicate(n = n, expr = bootstrap())
  
  mun_ests <- as.data.frame(sapply(1:n, FUN = function(x){unlist(reps[1,x])}))
  mun_ests$ID <- mun_list
  
  mun_levels <- as.data.frame(sapply(1:n, FUN = function(x){unlist(reps[2,x])}))
  mun_levels$ID <- mun_list
  
  mun_trust <-  as.data.frame(sapply(1:n, FUN = function(x){unlist(reps[3,x])}))
  mun_trust$ID <- mun_list
  
  comm_ests <- as.data.frame(sapply(1:n, FUN = function(x){unlist(reps[4,x])}))
  comm_ests$ID <- comm_list
  
  comm_levels <- as.data.frame(sapply(1:n, FUN = function(x){unlist(reps[5,x])}))
  comm_levels$ID <- comm_list
  
  comm_trust <- as.data.frame(sapply(1:n, FUN = function(x){unlist(reps[6,x])}))
  comm_trust$ID <- comm_list
  
  return(list(mun_ests, mun_levels, mun_trust, comm_ests, comm_levels, comm_trust))
}

# The following code takes about 5 minutes to run. The output is read in on line 123. Uncomment the code to conduct the bootstrapping.

# preds <- iterate(n =  5e2)

# mun_coef_ests <- preds[[1]]
# comm_coef_ests <- preds[[4]]
# mun_part_levs <- preds[[2]]
# comm_part_levs <- preds[[5]]
# mun_trust_levs <- preds[[3]]
# comm_trust_levs <- preds[[6]]


# mun_coef_ests <- left_join(mun_coef_ests, key)
# comm_coef_ests <- left_join(comm_coef_ests, key)
# mun_part_levels <- left_join(mun_part_levs, key)
# comm_part_levels <- left_join(comm_part_levs, key)
# mun_trust_levels <- left_join(mun_trust_levs, key)
# comm_trust_levels <- left_join(comm_trust_levs, key)

# save(mun_coef_ests, comm_coef_ests, mun_part_levels, comm_part_levels, 
#      mun_trust_levels, comm_trust_levels, file = "data/bootstrap_ests_lapop.Rdata")

load("data/bootstrap_ests_lapop.Rdata")

# Get estimates -- follows code analogous to bootstrap() but with realized sample

mun <- sapply(mun_list, FUN = function(x){
  summary(lm(mtg ~ mun_trust, data = lapop_sset, subset = ID == x))$coef[2,1]
})

mun_level <- tapply(lapop_sset$mtg, INDEX = lapop_sset$ID, FUN = mean, na.rm = T)
mun_level <- mun_level[names(mun_level) %in% mun_list]


comm <- sapply(comm_list, FUN = function(x){
  summary(lm(comm_assoc ~ comm_trust, data = lapop_sset, subset = ID == x))$coef[2,1]
})

comm_level <- tapply(lapop_sset$comm_assoc, INDEX = lapop_sset$ID, FUN = mean, na.rm = T)
comm_level <- comm_level[names(comm_level) %in% comm_list]

# col
mean(mun[grep(names(mun), pattern = "8_")]/mun_level[grep(names(mun), pattern = "8_")])
mean(comm[grep(names(comm), pattern = "8_")]/comm_level[grep(names(comm), pattern = "8_")])

# Select bootstrapped estimates from output and calculate 95% CIs
vars <- paste0("V",1:500)
# 95% CIs
mun_cis <- apply(sapply(vars, FUN = function(x){mun_coef_ests[,x]/mun_part_levels[,x]}), 1, quantile, probs = c(.025, .975))
comm_cis <- apply(sapply(vars, FUN = function(x){comm_coef_ests[,x]/comm_part_levels[,x]}), 1, quantile, probs = c(.025, .975))

# Figure 7 ----------------------------------------------------------------

lambda_df <- data.frame(lambda = c(mun/mun_level,comm/comm_level),
                       lo = c(mun_cis[1,], comm_cis[1,]),
                       hi = c(mun_cis[2,], comm_cis[2,]),
                       participation = c(rep("City Council", 115), rep("Community Association", 116)))%>%
  group_by(participation) %>%
  mutate(ranking = rank(lambda))

ggsave(lambda_df %>%
         mutate(sig = as.factor(1 * (lo > 0))) %>%
         ggplot(aes(x = ranking, y = lambda, color = sig)) + 
         facet_grid(~participation) +
         geom_hline(yintercept = 0, lty = 3, col = "dark red") + 
         scale_color_manual(values = c("dark grey", "black")) + 
         geom_point() + 
         geom_errorbar(aes(ymin = lo,
                           ymax = hi)) +
         scale_x_continuous("") + 
         ylab(expression(paste("Relationship between trust and attendance (", lambda, ")"))) + 
         theme_minimal() + 
         theme(legend.position = "n", axis.ticks.x = element_blank(),
               axis.text.x = element_blank()),
       file = "results/Figure_7.pdf", width = 9, height = 4)

# stats -- cited in paper
group_by(lambda_df, participation) %>%
  summarize(n = n(),
            pos_select = sum(lambda > 0),
            sig = sum(lo > 0 | hi < 0))


# Table A11 ---------------------------------------------------------------

# Calculate Pr(Attended | Trust = Highest) - Pr(Attended | Trust = Lowest) -- this is linearized
mun_actual <- sapply(mun_list, FUN = function(x){
  summary(lm(mtg ~ mun_trust, data = lapop_sset, subset = ID == x))$coef[2,1]
})
comm_actual <- sapply(mun_list, FUN = function(x){
  summary(lm(comm_assoc ~ comm_trust, data = lapop_sset, subset = ID == x))$coef[2,1]
})

# Calculate the probability of attendance Pr(Attended)
mun_lev_act <- sapply(mun_list, FUN = function(x){
  summary(lm(mtg ~1, data = lapop_sset, subset = ID == x))$coef[1,1]
})

comm_lev_act <- sapply(mun_list, FUN = function(x){
  summary(lm(comm_assoc ~ 1, data = lapop_sset, subset = ID == x))$coef[1,1]
})

# Calculate lambda: (Pr(Attended | Trust = Highest) - Pr(Attended | Trust = Lowest))/Pr(Attended)
mun_lambda_act <- mun_actual/mun_lev_act
comm_lambda_act <- comm_actual/comm_lev_act

# Calculate average level of trust in each institution
mun_trust_actual <- tapply(lapop_sset$mun_trust, INDEX = lapop_sset$ID, FUN = mean, na.rm = T)
mun_trust_actual <- mun_trust_actual[names(mun_trust_actual) %in% mun_list]

comm_trust_actual <- tapply(lapop_sset$comm_trust, INDEX = lapop_sset$ID, FUN = mean, na.rm = T)
comm_trust_actual <- comm_trust_actual[names(comm_trust_actual) %in% mun_list]

# Generates and outputs correlation matrix
print(xtable(round(cor(cbind(mun_actual, mun_lev_act, mun_lambda_act, mun_trust_actual, 
                       comm_actual, comm_lev_act, comm_lambda_act, comm_trust_actual)), 2)), 
      file = "results/Table_A11.tex")



# Figure A18 --------------------------------------------------------------
# calculate country averages
ctry_comm <- sapply(1:500, FUN = function(x) {
  lm(comm_coef_ests[,x]/comm_part_levels[,x] ~ -1+country, data = comm_coef_ests)$coef
})
ctry_mun <- sapply(1:500, FUN = function(x) {
  lm(mun_coef_ests[,x]/mun_part_levels[,x] ~ -1+country, data = mun_coef_ests)$coef
})

fig_a18a <- data.frame(country = sort(as.character(unique(key$pais))),
                 mean_theta = apply(ctry_comm, MAR = 1, mean),
                 ci_lo = apply(ctry_comm, MAR = 1, quantile, probs = .025),
                 ci_hi = apply(ctry_comm, MAR = 1, quantile, probs = .975),
                 measure = "Community Association") %>%
  ggplot(aes(x = mean_theta, y = reorder(country, mean_theta))) + 
  geom_point() + 
  facet_grid(~measure) + 
  geom_errorbarh(aes(xmin = ci_lo, xmax = ci_hi), height = 0) + 
  xlab(expression(paste("Relationship between trust and attendance (", gamma, ")"))) + 
  ylab("") + 
  theme_minimal() + 
  geom_vline(xintercept = 0, lty = 3, col = "red")


fig_a18b <- data.frame(country = sort(as.character(unique(key$pais))),
                 mean_theta = apply(ctry_mun, MAR = 1, mean),
                 ci_lo = apply(ctry_mun, MAR = 1, quantile, probs = .025),
                 ci_hi = apply(ctry_mun, MAR = 1, quantile, probs = .975),
                 measure = "City Council") %>%
  ggplot(aes(x = mean_theta, y = reorder(country, mean_theta))) + 
  geom_point() + 
  facet_grid(~measure) + 
  geom_errorbarh(aes(xmin = ci_lo, xmax = ci_hi), height = 0) + 
  xlab(expression(paste("Relationship between trust and attendance (", gamma, ")"))) + 
  ylab("") + 
  theme_minimal() + 
  geom_vline(xintercept = 0, lty = 3, col = "red")

p_patch <- fig_a18a + fig_a18b & xlab(NULL) & theme(plot.margin = margin(5.5, 5.5, 0, 5.5))

# Use the tag label as an x-axis label
fig_a18 <- wrap_elements(panel = p_patch) +
  labs(tag = expression(paste("Relationship between trust and attendance (", lambda, ")"))) +
  theme(
    plot.tag = element_text(size = rel(1)),
    plot.tag.position = "bottom"
  )

ggsave(fig_a18, 
       file = 'results/Figure_A18.pdf', width = 9, height =5)



